Functions
# Function to extract sample data as dataframe from phyloseq object
samdatAsDataframe <- function(ps) {
samdat <- phyloseq::sample_data(ps)
df <- data.frame(samdat, check.names = FALSE, stringsAsFactors = FALSE)
return(df)
}
# Function to rename variables in phyloseq object
ps_rename <- function(ps, ...) {
ps <- microViz::ps_get(ps)
df <- samdatAsDataframe(ps)
df <- dplyr::rename(.data = df, ...)
phyloseq::sample_data(ps) <- df
return(ps)
}
# SourceFolder function
source(here::here("Code", "R", "Functions", "StartFunctions", "sourceFolder.R"))
# Import all helper functions found in `/Functions`
sourceFolder(here::here("Code", "R", "Functions", "StartFunctions"), T)
##
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
##
## discard
## The following object is masked from 'package:readr':
##
## col_factor
## The following object is masked from 'package:coin':
##
## pvalue
## Loading required package: future
##
## Attaching package: 'future'
## The following object is masked from 'package:survival':
##
## cluster
##
## Attaching package: 'ggpubr'
## The following object is masked from 'package:microViz':
##
## stat_chull
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.18.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## ! pheatmap() has been masked by ComplexHeatmap::pheatmap(). Most of the arguments
## in the original pheatmap() are identically supported in the new function. You
## can still use the original function by explicitly calling pheatmap::pheatmap().
##
## Attaching package: 'ComplexHeatmap'
## The following object is masked from 'package:pheatmap':
##
## pheatmap
## Loading required package: data.table
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:lubridate':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year
## The following object is masked from 'package:purrr':
##
## transpose
## The following objects are masked from 'package:dplyr':
##
## between, first, last
## Loading required package: doParallel
## Loading required package: foreach
##
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
##
## accumulate, when
## Loading required package: iterators
## Loading required package: GUniFrac
## Registered S3 method overwritten by 'rmutil':
## method from
## print.response httr
## Loading required package: magrittr
##
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
##
## set_names
## The following object is masked from 'package:tidyr':
##
## extract
## Loading required package: rlang
##
## Attaching package: 'rlang'
## The following object is masked from 'package:magrittr':
##
## set_names
## The following object is masked from 'package:data.table':
##
## :=
## The following objects are masked from 'package:purrr':
##
## %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
## flatten_raw, invoke, splice
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4
##
## microbiome R package (microbiome.github.com)
##
##
##
## Copyright (C) 2011-2022 Leo Lahti,
## Sudarshan Shetty et al. <microbiome.github.io>
##
## Attaching package: 'microbiome'
## The following object is masked from 'package:vegan':
##
## diversity
## The following object is masked from 'package:scales':
##
## alpha
## The following object is masked from 'package:ggplot2':
##
## alpha
## The following object is masked from 'package:base':
##
## transform
## Loading required package: ape
##
## Attaching package: 'ape'
## The following object is masked from 'package:ggpubr':
##
## rotate
## The following object is masked from 'package:dplyr':
##
## where
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
##
## lmList
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## 3 files sourced from folder "/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Code/R/Functions/StartFunctions"
sourceFolder(here::here("Code", "R", "Functions", "HelperFunctions"), T)
## 9 files sourced from folder "/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Code/R/Functions/HelperFunctions"
# sourceFolder(here::here("Code", "R", "Functions", "AnalysisScripts"), T)
# Maaslin2
# Define a modified version of run_maaslin2
run_maaslin2_modified <- function(tmp.PS,
tmp.sample.names = "Sample",
tmp.rank = "Genus",
tmp.fixed,
tmp.reference,
tmp.output,
tmp.plot.scatter = T,
tmp.plot.heatmap = T,
tr = 'LOG',
norm = 'TSS'){
tmp.input.otu = tmp.PS %>% microViz::tax_agg(rank = "Genus") %>% otu_get() %>% as.data.frame() %>%
setNames(gsub(" ", "_", names(.)))
tmp.input.meta = microViz::samdat_tbl(tmp.PS) %>% as.data.frame()
row.names(tmp.input.meta) <- tmp.input.meta[[tmp.sample.names]]
# Create output directory if it doesn't exist
dir.create(tmp.output, recursive = TRUE, showWarnings = FALSE)
obj <- Maaslin2::Maaslin2(input_data = tmp.input.otu,
input_metadata = tmp.input.meta,
output = tmp.output,
min_prevalence = 0,
normalization = norm,
transform = tr,
fixed_effects = tmp.fixed,
reference = tmp.reference,
standardize = F,
plot_heatmap = tmp.plot.heatmap,
plot_scatter = tmp.plot.scatter,
cores = 8)
# Process results with corrected column selection
res <- obj$results %>%
filter(metadata == 'group' & value == 'case') %>%
mutate(df = N - length(tmp.fixed) - 1,
method = paste0('MaAsLin2 (', tr, ', ', norm, ')')) %>%
dplyr::select(feature, coef, stderr, df, pval, method) %>%
rename(taxon = feature,
est = coef,
se = stderr,
p = pval)
return(res)
}
#### GT TABLE STYLES
apply_GT_tabStyles <- function(gt_table, columns, coef_columns) {
for (i in seq_along(columns)) {
column <- columns[i]
coef_column <- coef_columns[i]
coef_vals <- gt_table$`_data`[[coef_column]]
colors <- gradient_palette(coef_vals)
for (j in seq_along(colors)) {
cell_color <- ifelse(is.na(coef_vals[j]), "grey95", colors[j])
cell_text <- "black" # ifelse(is.na(coef_vals[j]), "black", "white")
cell_size <- "medium"
bold_text <- ifelse(gt_table$`_data`[[column]][j] %in% c("+", "-"), "bold", "normal")
gt_table <- gt_table %>%
tab_style(
style = list(cell_fill(color = cell_color),
cell_text(color = cell_text,
weight = bold_text,
align = "center",
size = cell_size
)
),
locations = cells_body(
columns = !!sym(column),
rows = j
)
)
}
}
gt_table
}
# Function to create treatment-specific heatmap from MaAsLin2 results
create_treatment_heatmap <- function(maaslin_results, top_n = 50, title = "Treatment-specific Taxa Abundance Changes") {
# Set seed for reproducibility
set.seed(42)
# Filter for significant results and get top taxa
heatmap_data <- maaslin_results %>%
dplyr::filter(qval < 0.05) %>%
dplyr::arrange(qval) %>%
dplyr::slice_head(n = top_n) %>%
dplyr::select(Taxon, value, coef, Phylum) %>%
dplyr::distinct(Taxon, value, .keep_all = TRUE) # Remove duplicates if any
# Pivot to wide format with treatments as columns
heatmap_matrix <- heatmap_data %>%
tidyr::pivot_wider(
id_cols = c(Taxon, Phylum),
names_from = value,
values_from = coef,
values_fill = 0
) %>%
dplyr::arrange(Phylum, Taxon) %>%
tibble::column_to_rownames("Taxon")
# Extract phylum information for row annotations
phylum_info <- heatmap_matrix$Phylum
heatmap_matrix <- heatmap_matrix %>%
dplyr::select(-Phylum) %>%
as.matrix()
# Order columns according to treatment_order (only include treatments that exist in the data)
existing_treatments <- colnames(heatmap_matrix)
ordered_treatments <- treatment_order[treatment_order %in% existing_treatments]
heatmap_matrix <- heatmap_matrix[, ordered_treatments, drop = FALSE]
# Create row annotation data frame
annotation_row <- data.frame(
Phylum = phylum_info,
row.names = rownames(heatmap_matrix)
)
# Create color palette for phyla
phylum_colors <- setNames(
RColorBrewer::brewer.pal(length(unique(phylum_info)), "Set3"),
unique(phylum_info)
)
# Define annotation colors
ann_colors <- list(
Phylum = phylum_colors
)
# Create color palette for heatmap
# Use a diverging color palette centered at 0
max_abs_coef <- max(abs(heatmap_matrix), na.rm = TRUE)
heatmap_colors <- colorRampPalette(c("blue", "white", "red"))(100)
# Create the heatmap
pheatmap::pheatmap(
heatmap_matrix,
color = heatmap_colors,
breaks = seq(-max_abs_coef, max_abs_coef, length.out = 101),
annotation_row = annotation_row,
annotation_colors = ann_colors,
cluster_rows = TRUE,
cluster_cols = FALSE, # Keep treatment order as is
show_rownames = TRUE,
show_colnames = TRUE,
fontsize_row = 8,
fontsize_col = 10,
angle_col = 45,
main = title,
treeheight_row = 20,
treeheight_col = 0
)
}
# Function to create comparison tables for differentially abundant taxa
create_taxa_comparison_tables <- function(maaslin_results, section_title = "Taxa Comparison") {
# Set seed for reproducibility
set.seed(42)
if (is.null(maaslin_results) || nrow(maaslin_results) == 0) {
return(list(
summary_table = NULL,
overlap_table = NULL,
direction_table = NULL,
top10_table = NULL
))
}
# Filter for significant results only
sig_results <- maaslin_results %>%
dplyr::filter(qval < 0.05) %>%
dplyr::select(Taxon, value, coef, Phylum, qval)
if (nrow(sig_results) == 0) {
return(list(
summary_table = NULL,
overlap_table = NULL,
direction_table = NULL,
top10_table = NULL
))
}
# 1. Summary table - Number of significant taxa per treatment
summary_table <- sig_results %>%
dplyr::group_by(value) %>%
dplyr::summarise(
Total_Significant = dplyr::n(),
Positive_Effect = sum(coef > 0),
Negative_Effect = sum(coef < 0),
.groups = "drop"
) %>%
dplyr::arrange(factor(value, levels = treatment_order)) %>%
gt::gt() %>%
gt::tab_header(
title = paste("Summary of Significant Taxa by Treatment"),
subtitle = section_title
) %>%
gt::cols_label(
value = "Treatment",
Total_Significant = "Total Significant",
Positive_Effect = "Positive Effect",
Negative_Effect = "Negative Effect"
) %>%
gt::tab_style(
style = cell_text(weight = "bold"),
locations = cells_column_labels()
) %>%
gt::fmt_number(
columns = c(Total_Significant, Positive_Effect, Negative_Effect),
decimals = 0
)
# Color the Treatment column rows by treatment color
for (i in seq_along(treatment_order)) {
trt <- treatment_order[i]
if (trt %in% unique(sig_results$value)) {
trt_color <- treatment_color_scale[trt]
summary_table <- summary_table %>%
gt::tab_style(
style = cell_fill(color = trt_color),
locations = cells_body(
columns = "value",
rows = value == trt
)
) %>%
gt::tab_style(
style = cell_text(color = "white", weight = "bold"),
locations = cells_body(
columns = "value",
rows = value == trt
)
)
}
}
# 2. Top 10 taxa table - Most significant taxa by treatment
top10_data <- sig_results %>%
dplyr::group_by(value) %>%
dplyr::arrange(qval, .by_group = TRUE) %>%
dplyr::slice_head(n = 10) %>%
dplyr::ungroup() %>%
dplyr::select(Taxon, value, coef, Phylum, qval) %>%
dplyr::arrange(factor(value, levels = treatment_order), qval)
if (nrow(top10_data) > 0) {
top10_table <- top10_data %>%
gt::gt() %>%
gt::tab_header(
title = "Top 10 Most Significant Taxa by Treatment",
subtitle = paste(section_title, "- Ranked by q-value")
) %>%
gt::cols_label(
value = "Treatment",
coef = "Coefficient",
qval = "q-value"
) %>%
gt::tab_style(
style = cell_text(weight = "bold"),
locations = cells_column_labels()
) %>%
gt::fmt_scientific(
columns = qval,
decimals = 3
) %>%
gt::fmt_number(
columns = coef,
decimals = 3
) %>%
gt::tab_style(
style = cell_text(weight = "bold"),
locations = cells_body(columns = "Taxon")
)
# Color treatment rows by treatment color
for (i in seq_along(treatment_order)) {
trt <- treatment_order[i]
if (trt %in% unique(top10_data$value)) {
trt_color <- treatment_color_scale[trt]
top10_table <- top10_table %>%
gt::tab_style(
style = cell_fill(color = trt_color),
locations = cells_body(
columns = "value",
rows = value == trt
)
) %>%
gt::tab_style(
style = cell_text(color = "white", weight = "bold"),
locations = cells_body(
columns = "value",
rows = value == trt
)
)
}
}
} else {
top10_table <- NULL
}
# 3. Overlap table - Compare taxa between treatments
treatments <- unique(sig_results$value)
if (length(treatments) > 1) {
# Create list of taxa for each treatment
treatment_taxa <- list()
for (trt in treatments) {
treatment_taxa[[trt]] <- sig_results %>%
dplyr::filter(value == trt) %>%
dplyr::pull(Taxon)
}
# Create overlap matrix
overlap_matrix <- matrix(0, nrow = length(treatments), ncol = length(treatments))
rownames(overlap_matrix) <- treatments
colnames(overlap_matrix) <- treatments
for (i in seq_along(treatments)) {
for (j in seq_along(treatments)) {
if (i == j) {
overlap_matrix[i, j] <- length(treatment_taxa[[treatments[i]]])
} else {
overlap_matrix[i, j] <- length(intersect(treatment_taxa[[treatments[i]]],
treatment_taxa[[treatments[j]]]))
}
}
}
# Convert to data frame for GT table
overlap_df <- as.data.frame(overlap_matrix) %>%
tibble::rownames_to_column("Treatment") %>%
dplyr::arrange(factor(Treatment, levels = treatment_order))
overlap_table <- overlap_df %>%
gt::gt() %>%
gt::tab_header(
title = "Taxa Overlap Between Treatments",
subtitle = paste(section_title, "- Number of shared significant taxa")
) %>%
gt::tab_style(
style = cell_text(weight = "bold"),
locations = cells_column_labels()
) %>%
gt::fmt_number(
columns = -Treatment,
decimals = 0
)
# Apply treatment colors to column headers and row headers
for (i in seq_along(treatments)) {
trt <- treatments[i]
if (trt %in% names(treatment_color_scale)) {
trt_color <- treatment_color_scale[trt]
# Color column header
overlap_table <- overlap_table %>%
gt::tab_style(
style = cell_fill(color = trt_color),
locations = cells_column_labels(columns = trt)
) %>%
gt::tab_style(
style = cell_text(color = "white", weight = "bold"),
locations = cells_column_labels(columns = trt)
)
# Color row header (first column) for matching treatment
overlap_table <- overlap_table %>%
gt::tab_style(
style = cell_fill(color = trt_color),
locations = cells_body(
columns = "Treatment",
rows = Treatment == trt
)
) %>%
gt::tab_style(
style = cell_text(color = "white", weight = "bold"),
locations = cells_body(
columns = "Treatment",
rows = Treatment == trt
)
)
}
}
} else {
overlap_table <- NULL
}
# 4. Direction comparison table - Compare positive/negative effects
if (length(treatments) > 1) {
direction_data <- sig_results %>%
dplyr::group_by(value) %>%
dplyr::summarise(
Positive_Taxa = list(Taxon[coef > 0]),
Negative_Taxa = list(Taxon[coef < 0]),
.groups = "drop"
)
# Create direction comparison matrix
direction_matrix <- matrix("", nrow = length(treatments), ncol = length(treatments))
rownames(direction_matrix) <- treatments
colnames(direction_matrix) <- treatments
for (i in seq_along(treatments)) {
for (j in seq_along(treatments)) {
if (i == j) {
# For diagonal cells, show total positive and negative counts for that treatment
# Find the correct treatment data
trt_name <- treatments[i]
trt_data <- direction_data[direction_data$value == trt_name, ]
if (nrow(trt_data) > 0) {
pos_count <- length(trt_data$Positive_Taxa[[1]])
neg_count <- length(trt_data$Negative_Taxa[[1]])
direction_matrix[i, j] <- paste0(
"Pos: ", pos_count,
" | Neg: ", neg_count
)
} else {
direction_matrix[i, j] <- "Pos: 0 | Neg: 0"
}
} else {
# Count taxa with same direction between two different treatments
pos_i <- direction_data$Positive_Taxa[[i]]
neg_i <- direction_data$Negative_Taxa[[i]]
pos_j <- direction_data$Positive_Taxa[[j]]
neg_j <- direction_data$Negative_Taxa[[j]]
# Find taxa that are positive in both treatments
same_pos <- length(intersect(pos_i, pos_j))
# Find taxa that are negative in both treatments
same_neg <- length(intersect(neg_i, neg_j))
direction_matrix[i, j] <- paste0(
"Same Pos: ", same_pos,
" | Same Neg: ", same_neg
)
}
}
}
direction_df <- as.data.frame(direction_matrix) %>%
tibble::rownames_to_column("Treatment") %>%
dplyr::arrange(factor(Treatment, levels = treatment_order))
direction_table <- direction_df %>%
gt::gt() %>%
gt::tab_header(
title = "Direction of Effect Comparison",
subtitle = paste(section_title, "- Shared taxa with same direction of effect")
) %>%
gt::tab_style(
style = cell_text(weight = "bold"),
locations = cells_column_labels()
)
# Apply treatment colors to column headers and row headers
for (i in seq_along(treatments)) {
trt <- treatments[i]
if (trt %in% names(treatment_color_scale)) {
trt_color <- treatment_color_scale[trt]
# Color column header
direction_table <- direction_table %>%
gt::tab_style(
style = cell_fill(color = trt_color),
locations = cells_column_labels(columns = trt)
) %>%
gt::tab_style(
style = cell_text(color = "white", weight = "bold"),
locations = cells_column_labels(columns = trt)
)
# Color row header (first column) for matching treatment
direction_table <- direction_table %>%
gt::tab_style(
style = cell_fill(color = trt_color),
locations = cells_body(
columns = "Treatment",
rows = Treatment == trt
)
) %>%
gt::tab_style(
style = cell_text(color = "white", weight = "bold"),
locations = cells_body(
columns = "Treatment",
rows = Treatment == trt
)
)
}
}
} else {
direction_table <- NULL
}
return(list(
summary_table = summary_table,
top10_table = top10_table,
overlap_table = overlap_table,
direction_table = direction_table
))
}
Import Data
ps.tmp <- readRDS("/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Data/Robjects/pseq_uncleaned_05052025.rds")
ps.cleaned <-
ps.tmp %>%
## Update Metadata
ps_rename(Time = Timepoint) %>%
microViz::ps_mutate(
Treatment = case_when(
Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "A- T- P-",
Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "A- T- P+",
Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "A+ T- P-",
Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "A+ T- P+",
Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "A- T+ P-",
Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "A- T+ P+",
Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "A+ T+ P-",
Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "A+ T+ P+",
TRUE ~ "Unknown"
), .after = "Pathogen"
) %>%
microViz::ps_mutate(Sample = fecal.sample.number, .before = 1) %>%
microViz::ps_mutate(Sample = gsub("^f", "", Sample)) %>%
microViz::ps_filter(Treatment != "Unknown") %>%
microViz::ps_mutate(
History = case_when(
Antibiotics + Temperature == 0 ~ 0,
Antibiotics + Temperature == 1 ~ 1,
Antibiotics + Temperature == 2 ~ 2,
), .after = "Treatment"
) %>%
## Additional metadata updates, factorizing metadata
microViz::ps_mutate(
# Create treatment code
treatment_code = case_when(
Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "Aneg_Tneg_Pneg",
Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Aneg_Tneg_Ppos",
Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Apos_Tneg_Pneg",
Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Apos_Tneg_Ppos",
Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Aneg_Tpos_Pneg",
Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Aneg_Tpos_Ppos",
Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Apos_Tpos_Pneg",
Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Apos_Tpos_Ppos"
),
# Create treatment group factor
treatment_group = case_when(
Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Parasite",
Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Antibiotics",
Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Antibiotics_Parasite",
Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Temperature",
Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Temperature_Parasite",
Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Antibiotics_Temperature",
Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Antibiotics_Temperature_Parasite",
TRUE ~ "Control"
),
# Convert to factor with appropriate levels
treatment_group = factor(treatment_group,
levels = c("Control", "Parasite",
"Antibiotics", "Antibiotics_Parasite",
"Temperature", "Temperature_Parasite",
"Antibiotics_Temperature", "Antibiotics_Temperature_Parasite")
),
treatment_code = factor(treatment_code, levels = treatment_order),
# Create time point factor
time_point = factor(Time, levels = c(0, 14, 18, 25, 29, 60)),
# Create pathogen status factor
pathogen_status = factor(ifelse(Pathogen == 1, "Exposed", "Unexposed"),
levels = c("Unexposed", "Exposed")),
# Create sex factor
sex = factor(Sex, levels = c("M", "F"))
) %>%
microViz::ps_mutate(Treatment = factor(Treatment, levels = treatment_order)) %>%
microViz::ps_mutate(Exp_Type = case_when(
Treatment %in% c("A- T- P-", "A- T- P+") ~ "No prior stressor(s)",
Treatment %in% c("A+ T- P-", "A+ T- P+") ~ "Antibiotics",
Treatment %in% c("A- T+ P-", "A- T+ P+") ~ "Temperature",
Treatment %in% c("A+ T+ P-", "A+ T+ P+") ~ "Combined",
)) %>%
microViz::ps_mutate(Exp_Type = factor(Exp_Type, levels = c("No prior stressor(s)", "Antibiotics", "Temperature", "Combined"))) %>%
# Fix names for taxonomic ranks not identified
microViz::tax_fix(suffix_rank = "current", anon_unique = T, unknown = NA) %>%
# Filter for any samples that contain more than 5000 reads
microViz::ps_filter(sample_sums(.) > 5000) %>%
# Any taxa not found in at least 3 samples are removed
microViz::tax_filter(min_prevalence = 3, undetected = 0) %>%
# Remove any unwanted reads
microViz::tax_select(c("Mitochondria", "Chloroplast", "Eukaryota"), deselect = TRUE) %>%
microViz::tax_select(c("Bacteria, Phylum"), deselect = TRUE)
# ps.cleaned %>% microViz::samdat_tbl()
